function [mom, panel_out] = compute_moments(eq_SS, param, glob, options)

    %% markups
    eq_SS.mc     = param.omega*1./exp(glob.sf(:,2));
    eq_SS.markup = eq_SS.pf./eq_SS.mc;
    eq_SS.sales  = eq_SS.pf.*eq_SS.yf;
        
    eq_SS.costs  = eq_SS.l*param.omega;

    
    mom.mu_sw    = sum(eq_SS.sales.*eq_SS.L.*eq_SS.markup)/sum(eq_SS.sales.*eq_SS.L);
    mom.mu_cw    = sum(eq_SS.costs.*eq_SS.L.*eq_SS.markup)/sum(eq_SS.costs.*eq_SS.L);

    mom.mu_uw    = sum(eq_SS.L.*eq_SS.markup)/sum(eq_SS.L);
    
    q            = eq_SS.yf;
    elast        = param.sigma*q.^(-param.epsilon/param.sigma);
    mu_sigma     = elast./(elast-1);
    mom.cw_sigma = sum(mu_sigma.*eq_SS.costs.*eq_SS.L)/sum(eq_SS.costs.*eq_SS.L);
    
    mom.mu_disp  = sum((eq_SS.markup - mom.mu_uw).^2.*eq_SS.L);
    
    mom.mu_sales = sum(eq_SS.sales.*eq_SS.L);
    
    mom.sales_disp  = sum((eq_SS.sales - mom.mu_sales).^2.*eq_SS.L);
    
    param.w = 1; param.D = eq_SS.D; param.C = 1;
    
    acr             =  menufun('C',glob.sf,eq_SS.yf,param,glob,options); 
    
    mom.cost_size = sum(eq_SS.L.*acr)/sum(eq_SS.L.*eq_SS.sales);
    
    %%
     Lf = reshape(eq_SS.L, glob.Nsf,1);
     F = randsample(size(Lf,1), 1e6, true, Lf);

     prices  = eq_SS.pf(F);
     sales   = eq_SS.yf(F).*prices;
          
    prices  = eq_SS.pf(F);
    sales   = eq_SS.yf(F).*prices;

     
     prods   = glob.sf(F,2);

     markups = prices.*exp(prods);
     costs   = sales./prices./exp(prods);
 

    %% top 1% and 5% sales share
    mom.top10 = sum(maxk(sales, floor(.10*length(sales))))/sum(sales);
    mom.top05 = sum(maxk(sales, floor(.05*length(sales))))/sum(sales);
    mom.top01 = sum(maxk(sales, floor(.01*length(sales))))/sum(sales);

    mom.top01_emp = sum(maxk(costs, floor(.01*length(sales))))/sum(costs);

    mom.bottom87 = sum(mink(sales, floor(.877*length(sales))))/sum(sales); %15.4 percent of sales
    mom.frac_rel_sales_below_1 = sum(sales < mean(sales))/length(sales);
    mom.frac_rel_sales_below_10 = sum(sales < 10*mean(sales))/length(sales);
    mom.frac_rel_sales_below_50 = sum(sales < 50*mean(sales))/length(sales);

     %% simulate panel
    panel = ddpsimul(eq_SS.Q, F(1:100000),1);
  
%%
    
    prices  = eq_SS.pf(panel);

    prods   = glob.sf(:,2);
    prods   = prods(panel);
    labor   = eq_SS.yf(panel)./exp(prods);
    sales   = eq_SS.yf(panel).*prices;
        
    % I want to pick a sample that "looks like" Compustat.
    % If I pick big firms in first period, in 2nd period they're too small
    % because of mean reversion.
    % Instead, pick top firms among sum of total output in 1st and second
    % period
    

    sales   = sum(sales,2);
    labor   = sum(labor,2);
    [~,y] = sort(sales, 'descend');
    
    sample = labor(y);
    idx    = floor(length(sample)*.01); find(cumsum(sample) >= .3*sum(sample), 1);
    idx    = y(1:idx);
    
    sales   = eq_SS.yf(panel).*prices;
    labor   = eq_SS.yf(panel)./exp(prods);

    panel = panel(idx,:);
    prices  = prices(idx,:);
    prods   = prods(idx,:);
    sales   = sales(idx,:);
    markups = prices.*exp(prods);
    
    labor = labor(idx,:);
    
    mom.labor_churn = var(diff(log(labor),[],2));
    mom.sales_churn = var(diff(log(sales),[],2));
    mom.markups_churn = var(diff(log(markups),[],2));
    
    panel_out = [];
    panel_out.panel = panel;
    panel_out.labor = labor;
    panel_out.sales = sales;
    panel_out.markups = markups;
    panel_out.prices = prices;
    panel_out.prods = prods;
   
    
    %% run markups regression within firm
    %  also do labor regression
    
    %% compute slope coefficient
    lm = fitlm(log(sales(:,2)) - log(sales(:,1)),log(markups(:,2)) - log(markups(:,1)));
    mom.beta_mu = lm.Coefficients{2,1};
    
    lm = fitlm(log(sales(:,2)) - log(sales(:,1)),log(labor(:,2)) - log(labor(:,1)));
    mom.beta_l = lm.Coefficients{2,1};

    gs = (sales(:,2) - sales(:,1))./(.5*sales(:,2) + .5*sales(:,1));
    gl = (labor(:,2) - labor(:,1))./(.5*labor(:,2) + .5*labor(:,1));

    lm = fitlm(gs, gl);
    mom.beta_l_2 = lm.Coefficients{2,1};


    %%
    lm = fitlm(log(sales(:)),log(labor(:)));
    mom.beta_l_btw = lm.Coefficients{2,1};


    %% labor adjustment among large firms
    
    mom.large_jcjd10 = sum(diff(log(labor),[],2) > .1)/size(labor,1); % this moment is 28% in compustat (sample: post 2010)
    
    %% job creation/job destruction rates
    
    dL = (eq_SS.l - glob.sf(:,1))./glob.sf(:,1);
    %dL(glob.sf(:,1) == 0) = 0;
        
    mom.jc10 = sum(eq_SS.L.*(dL < .1).*(dL > 0))/sum(eq_SS.L);  % target: ~.05 quarterly, so ~.0125?
    mom.jd10 = sum(eq_SS.L.*(dL > -.1).*(dL < 0))/sum(eq_SS.L); % target: ~.05
    
    %% relative employment of exiters
    
    av_size_exiter = (sum(param.gamma*eq_SS.l.*eq_SS.L) + sum(eq_SS.l.*eq_SS.L.*eq_SS.p_exit))/sum(param.gamma*eq_SS.L+ eq_SS.L.*eq_SS.p_exit);
    
    av_labor = sum(eq_SS.l.*eq_SS.L)/sum(eq_SS.L);
    
    mom.rel_size_exiter = av_size_exiter/av_labor;
    
    %% Relative size of entrants vs. incumbents
    
    mom.entry_rate = eq_SS.mass/sum(eq_SS.L);    
    
    G = eq_SS.G;
    G =  param.M*eq_SS.Q'*G;

    av_size_entrant = sum(eq_SS.l.*G.*(1-eq_SS.p_exit))/sum(G.*(1-eq_SS.p_exit));
    
    Lincumbent = eq_SS.L - eq_SS.mass*G;
    Lincumbent = Lincumbent.*(1-eq_SS.p_exit);
    av_size_incumbent = sum(eq_SS.l.*Lincumbent)/sum(Lincumbent);
    
    mom.size_ratio = av_size_entrant/av_size_incumbent;
    mom.av_size_incumbent = av_size_incumbent;
    
    %% Share of employment among young firms
    
    % Entrants: 5.5% in 2006, 4% in 2014
    % 5 and younger:  30% in 2006, 22% in 2014
    
    mom.labor_demand      = sum(eq_SS.l.*eq_SS.L);
    
    mom.entrant_emp_share = sum(param.M.*eq_SS.l.*eq_SS.G)./mom.labor_demand;
    
    % get distribution of firms of ages 1, 2, 3, 4, 5
     a1      = param.M*eq_SS.Q'*((1-param.gamma)*(1-eq_SS.p_exit).*eq_SS.G);
     a2      = eq_SS.Q'*((1-param.gamma)*(1-eq_SS.p_exit).*a1);
     a3      = eq_SS.Q'*((1-param.gamma)*(1-eq_SS.p_exit).*a2);
     a4      = eq_SS.Q'*((1-param.gamma)*(1-eq_SS.p_exit).*a3);
     a5      = eq_SS.Q'*((1-param.gamma)*(1-eq_SS.p_exit).*a4);
     a6      = eq_SS.Q'*((1-param.gamma)*(1-eq_SS.p_exit).*a5);
% 
    mom.employment_alt5 = sum(eq_SS.l.*G) + ...
                          sum(eq_SS.l.*a1) + ...
                          sum(eq_SS.l.*a2) + ...
                          sum(eq_SS.l.*a3) + ...
                          sum(eq_SS.l.*a4) + ...
                          sum(eq_SS.l.*a5);

    mom.young_emp_share = mom.employment_alt5/mom.labor_demand;
    
    rev = eq_SS.pf.*eq_SS.yf;

    
    mom.revenue_alt5    = sum(rev.*G) + ...
                          sum(rev.*a1) + ...
                          sum(rev.*a2) + ...
                          sum(rev.*a3) + ...
                          sum(rev.*a4) + ...
                          sum(rev.*a5);
    
    mom.revenue_a1     = sum(rev.*a1);                  
                      
    %% Life cycle of the firm
%     
    mom.emp_a0 = sum(eq_SS.l.*G)/sum(G);
    mom.emp_a1 = sum(eq_SS.l.*a1)/sum(a1);
    mom.emp_a2 = sum(eq_SS.l.*a2)/sum(a2);
    mom.emp_a3 = sum(eq_SS.l.*a3)/sum(a3);
    mom.emp_a4 = sum(eq_SS.l.*a4)/sum(a4);
    mom.emp_a5 = sum(eq_SS.l.*a5)/sum(a5);
    mom.emp_a6 = sum(eq_SS.l.*a6)/sum(a6);
    
    mom.emp_av = sum(eq_SS.l.*eq_SS.L)/sum(eq_SS.L);
    
    mom.prod_a0 = sum(glob.sf(:,2).*G)/sum(G);
    mom.prod_a1 = sum(glob.sf(:,2).*a1)/sum(a1);
    mom.prod_a2 = sum(glob.sf(:,2).*a2)/sum(a2);
    mom.prod_a3 = sum(glob.sf(:,2).*a3)/sum(a3);
    mom.prod_a4 = sum(glob.sf(:,2).*a4)/sum(a4);
    mom.prod_a5 = sum(glob.sf(:,2).*a5)/sum(a5);
    mom.prod_a6 = sum(glob.sf(:,2).*a6)/sum(a6);
    
    mom.prod_av = sum(glob.sf(:,2).*eq_SS.L)/sum(eq_SS.L);

    
    rev = eq_SS.pf.*eq_SS.yf;
    
    mom.rev_a0 = sum(rev.*G)/sum(G);
    mom.rev_a1 = sum(rev.*a1)/sum(a1);
    mom.rev_a2 = sum(rev.*a2)/sum(a2);
    mom.rev_a3 = sum(rev.*a3)/sum(a3);
    mom.rev_a4 = sum(rev.*a4)/sum(a4);
    mom.rev_a5 = sum(rev.*a5)/sum(a5);
    mom.rev_a6 = sum(rev.*a6)/sum(a6);
    
    mom.rev_av = sum(rev.*eq_SS.L)/sum(eq_SS.L);

    mkp = eq_SS.pf.*exp(glob.sf(:,2));
    
    mom.mkp_a0 = sum(mkp.*G.*eq_SS.l)/sum(G.*eq_SS.l);
    mom.mkp_a1 = sum(mkp.*a1.*eq_SS.l)/sum(a1.*eq_SS.l);
    mom.mkp_a2 = sum(mkp.*a2.*eq_SS.l)/sum(a2.*eq_SS.l);
    mom.mkp_a3 = sum(mkp.*a3.*eq_SS.l)/sum(a3.*eq_SS.l);
    mom.mkp_a4 = sum(mkp.*a4.*eq_SS.l)/sum(a4.*eq_SS.l);
    mom.mkp_a5 = sum(mkp.*a5.*eq_SS.l)/sum(a5.*eq_SS.l);
    mom.mkp_a6 = sum(mkp.*a6.*eq_SS.l)/sum(a6.*eq_SS.l);

    mom.mkp_av = sum(mkp.*eq_SS.L.*eq_SS.l)/sum(eq_SS.L.*eq_SS.l);

    %% distribution of relative sales
    % whole economy statistics
         Lf = reshape(eq_SS.L, glob.Nsf,1);
     F = randsample(size(Lf,1), 10000, true, Lf);

     
    panel = ddpsimul(eq_SS.Q, F,2);
    
    labor = eq_SS.l(panel);
    
    sales = eq_SS.p(panel).*eq_SS.y(panel);
    
    d_labor = [labor(:,2) - labor(:,1), labor(:,3) - labor(:,2)];
    d_labor = d_labor./(.5*labor(:,1:2) + .5*labor(:,2:3));
    
    d_sales = [sales(:,2) - sales(:,1), sales(:,3) - sales(:,2)];
    d_sales = d_sales./(.5*sales(:,1:2) + .5*sales(:,2:3));
 
    
    mom.reallocation_rate = mean(mean(abs(d_labor)));
    
    labor_moments   = cov(d_labor); 
   	sales_moments   = cov(d_sales); 
    
    mom.demp_auto   = labor_moments(1,2)/labor_moments(1,1);
    mom.demp_var    = labor_moments(1,1);
    
    mom.dsales_auto = sales_moments(1,2)/labor_moments(1,1);
    mom.dsales_var  = sales_moments(1,1);
    
 %   mom.hiring_rate       = mean(d_labor);
    
    
    
end

